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ABSTRACT 


The  report  describes  the  status  and  results  at  the  end  of  the  second  year  of  a  project  to  reduce  the 
magnitude  threshold  for  which  surface  waves  can  be  identified  and  measured  reliably,  and  to 
improve  the  accuracy  of  surface  wave  measurement,  using  phase-matched  filtering,  development 
of  global  regionalized  earth  and  dispersion  models,  and  other  techniques.  In  year  two,  we  have 
focused  on  two  topics:  improvements  to  global  earth  models  and  dispersion  maps,  and  improved 
techniques  for  measuring  surface  wave  amplitudes.  We  also  completed  work  on  implementation 
and  testing  of  azimuth  estimation  techniques  at  three  component  stations  based  on  polarization 
analysis.  Global  regionalized  earth  and  dispersion  models  are  being  developed  by  inversion  of  a 
very  large  data  set  of  phase  and  group  velocity  dispersion  measurements.  The  complete  data  set 
now  contains  over  one  million  dispersion  data  points.  The  dispersion  measurements  are  inverted 
for  earth  structure,  and  the  earth  structure  is  then  used  to  generate  dispersion  predictions.  A 
significant  improvement  in  the  inversion  procedure  over  the  past  year  has  been  the  introduction 
of  the  capability  to  vary  damping  and  smoothing  parameters  for  each  model.  This  allows  us  to 
improve  the  data  fit  for  many  models,  while  still  retaining  realistic  earth  models  in  all  areas.  In 
addition  to  improving  earth  and  dispersion  models,  we  have  implemented  and  tested  procedures 
for  measuring  surface  wave  amplitudes  at  short  periods  (5-15  seconds)  and  at  regional  distances. 
We  find  that  the  best  procedure  is  to  measure  surface  wave  amplitudes  at  periods  greater  than  12 
seconds  at  all  distance  ranges  for  three  reasons:  1)  earthquake  spectra  tend  to  decrease  with 
frequency,  degrading  discrimination;  2)  unlike  time  domain  measurements,  spectral 
measurements  can  be  made  accurately  at  lower  frequencies  at  close  distances;  and  3)  S/N 
remains  higher  at  periods  greater  than  12  seconds  even  for  very  close  distances. 
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EXECUTIVE  SUMMARY 


Surface  waves  are  of  primary  importance  for  nuclear  monitoring  because  the  Ms:mb  discriminant 
and  its  regional  variants  are  among  the  most  reliable  means  of  determining  whether  an  event  is 
an  earthquake  or  an  explosion.  The  primary  goal  of  this  project  is  to  reduce  the  magnitude 
threshold  for  which  surface  waves  can  be  identified  and  measured  reliably  and  to  improve  the 
accuracy  of  surface  wave  measurement  using  phase-matched  filtering  and  global  regionalized 
earth  and  dispersion  models. 

Global  regionalized  earth  and  dispersion  models  are  being  developed  by  inversion  of  a  very  large 
data  set  of  phase  and  group  velocity  dispersion  measurements.  The  complete  data  set  now 
contains  over  one  million  dispersion  data  points.  The  dispersion  measurements  are  inverted  for 
earth  structure,  and  the  earth  structure  is  then  used  to  generate  dispersion  predictions.  This  is 
accomplished  in  the  following  way.  The  inversion  is  performed  for  approximately  600  distinct 
base  structures,  which  were  originally  derived  from  the  Crust  2.0  models  over  the  AK135  mantle 
model.  The  Moho  depth,  bathymetry,  and  sediment  depths  vary  on  a  one-degree  grid.  Moho 
depths  are  derived  from  Crust  2.0,  sediment  depths  from  the  Laske  and  Masters  sediment  maps, 
and  bathymetry  from  ETOP05.  Moho  depth,  bathymetry,  and  sediment  properties  are  fixed  in 
the  inversion,  while  crust  and  upper  mantle  velocities  are  allowed  to  vary  in  the  base  models. 
Phase  and  group  velocity  dispersion  curves  are  calculated  for  each  of  64,800  models  on  the  one- 
degree  grid.  The  phase  velocity  dispersion  curves  are  then  used  to  calculate  phase-matched  filters 
to  improve  detection. 

One  of  the  difficulties  of  performing  such  a  large,  heterogeneous  inversion  is  finding  optimum 
values  for  smoothing  and  damping  (regularizing)  parameters.  This  is  important  because  too 
much  smoothing/damping  will  increase  data  misfit,  and  too  little  will  produce  unrealistic  earth 
models.  In  this  case,  smoothing  minimizes  the  gradient  of  each  structure  between  specified 
discontinuities,  while  damping  minimizes  the  difference  between  the  model  and  the  starting 
model.  Discontinuities  occur  at  the  Moho  and  at  the  base  of  the  surface  sediments.  In  a  few  cases 
where  there  is  sufficient  high  frequency  information  to  resolve  shallower  structure,  inversion  is 
performed  for  deeper  sediments,  which  introduces  another  discontinuity.  A  significant 
improvement  in  the  inversion  procedure  over  the  past  year  has  been  the  introduction  of  the 
capability  to  vary  damping  and  smoothing  parameters  for  each  model.  This  allows  us  to  improve 
the  data  fit  for  many  models,  while  still  retaining  realistic  earth  models  in  all  areas.  We  are 
performing  2D  inversions  (inversion  of  dispersion  at  discrete  frequencies  to  form  spatial 
dispersion  maps)  to  identify  regions  where  additional  parameterization  is  needed  in  the  3D 
inversions.  The  additional  parameterization  takes  the  form  of  introduction  of  new  base  structures, 
merging  of  base  structures,  or  adjustments  to  boundaries  between  base  structures. 

In  addition  to  improving  earth  and  dispersion  models,  we  have  implemented  and  tested 
procedures  for  measuring  surface  wave  amplitudes  at  short  periods  (5-15  seconds)  and  at 
regional  distances.  We  are  identifying  optimum  procedures  for  measuring  path  corrected  spectral 
magnitudes  and  then  comparing  the  results  with  other  procedures,  such  as  Marshall-Basham 
amplitude  corrections,  which  were  designed  to  correct  time  domain  amplitudes  for  frequency  and 
structure  dependence.  Path  corrected  spectral  magnitudes  should  be  independent  of  distance  and 
only  weakly  dependent  on  frequency  for  shallow  events.  At  higher  frequencies  and  longer 
distances  the  amplitude  correction  depends  on  having  accurate  Q  models.  We  find  that  the  best 
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procedure  is  to  measure  surface  wave  amplitudes  at  periods  greater  than  12  seconds  at  all 
distance  ranges  for  three  reasons:  1)  earthquake  spectra  tend  to  decrease  with  frequency, 
degrading  discrimination;  2)  unlike  time  domain  measurements,  spectral  measurements  can  be 
made  accurately  at  lower  frequencies  at  close  distances;  and  3)  S/N  remains  higher  at  periods 
greater  than  12  seconds  even  for  very  close  distances. 

We  evaluate  the  performance  of  a  new  algorithm  for  determining  Rayleigh  wave  propagation 
direction  using  the  Rayleigh  wave  polarization.  The  algorithm  uses  the  correlation  of  the 
horizontal  and  Hilbert  transformed  vertical  seismograms  to  estimate  backazimuth,  where  the 
horizontal  component  is  rotated  through  the  full  range  of  possible  backazimuths.  The  vertical  and 
radial  seismograms’  correlation  coefficient  predicts  the  accuracy  of  the  backazimuth  estimate. 
The  effect  on  backazimuth  residuals  of  group  velocity  window,  passband,  and  different  measures 
of  peak  coirelation,  is  evaluated  using  a  large  data  set.  We  find  that  the  algorithm  provides 
accurate  estimates  of  Rayleigh  wave  polarization  and  is  a  significant  improvement  over  current 
automatic  processing  procedures.  The  algorithm  can  be  combined  with  a  dispersion  test  to 
improve  detection  capability. 
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2. .  GLOBAL  EARTH  MODELS  AND  SURFACE  WAVE  DISPERSION  MAPS 


To  improve  the  detection  and  measurement  of  surface  waves  it  is  important  to  make  good 
dispersion  predictions.  Dispersion  is  used  in  two  ways:  1)  to  construct  phase-matched  filters  to 
improve  signal  to  noise  ratio,  and  2)  to  provide  sufficient  time  resolution  to  test  that  the  surface 
wave  is  correctly  associated  with  a  particular  event.  In  our  work,  surface  wave  dispersion 
predictions  are  based  on  dispersion  measurements  for  ray  paths  from  all  over  the  globe.  We 
make  these  predictions  via  models  of  the  earth.  Alternate  methods  could  be  developed  which 
would  depend  on  interpolation  schemes,  such  as  kriging  and  which  would  not  make  use  of  earth 
models.  An  important  advantage  gained  from  using  earth  models  is  that  we  can  include 
information  from  other  studies  leading  to  physically  reasonable  constraints  on  dispersion.  For 
our  earth  models  this  information  consists  of  the  boundaries  between  geologic  zones,  bathymetry 
of  oceans,  thicknesses  of  sediments  and  ice,  Moho  depths,  and  prior  estimates  of  seismic 
velocities  derived  from  Crust  2.0  and  AK  135  earth  models.  These  constraints  are  especially 
important  for  filling  in  the  gaps  found  in  the  path  coverage  of  our  data  set  and  they  enable 
prediction  along  paths  unlike  any  of  the  paths  in  the  data  set.  We  perform  non-linear  least 
squares  inversion  of  the  dispersion  data  for  two  types  of  models.  First  are  3D  earth  models  where 
the  adjustable  parameters  are  the  shear  wave  velocities  of  layers.  Second  are  2D  group  velocity 
models  determined  by  inverting  dispersion  measurements  made  in  narrow  frequency  bands.  The 
2D  models  are  used  as  a  guide  in  the  parameterization  of  the  3D  earth  model. 

2.1  Description  of  3D  earth  model 

The  3D  earth  model  is  described  briefly  here  and  in  greater  detail  in  Stevens  et  al.  (2002).  It 
consists  of  l°xl°  blocks  and  is  made  up  of  layers  of  ice,  water,  sediments,  crust  and  upper 
mantle.  Currently  this  model  depends  on  8918  free  parameters  which  are  adjusted  by  least 
squares  fitting  to  Rayleigh  wave  dispersion  data.  The  free  parameters  are  the  S-wave  velocities 
of  layers  of  572  different  model  types.  Other  constrained  parameters  in  the  model  are  P  wave 
velocities,  densities  and  Q.  The  model  types  are  based  on  the  Crust  2.0  2°x2°  crustal  tjq)®^ 
(Bassin  et  al.,  2000  and  Laske  et  al.  2001)  and  also  on  ocean  ages  (Stevens  and  Adams,  2000). 
The  top  few  km  of  the  model  (consisting  of  water,  ice  and/or  sediments)  are  fixed  and  match  data 
from  one  degree  bathymetry  maps  made  by  averaging  Etopo5  five  minute  measurements  of 
topography,  and  Laske  and  Masters  (1997)  1  degree  maps  of  sediments.  There  is  an  explicit 
discontinuity  between  the  bottom  of  the  sediments  and  the  crust.  There  are  three  or  more  crustal 
layers.  The  Crust  2.0  models  which  were  the  starting  point  for  these  structures  have  three  cmstal 
layers,  but  we  found  it  necessary  to  add  more  layers  in  regions  of  thick  crust.  There  is  another 
explicit  discontinuity  at  the  Crust/Mantle  boundary.  The  Moho  depth  is  derived  from  Crust  2.0 
and  varies  on  a  2“  grid.  The  mantle  starting  model  is  derived  from  AK135  (Kennett,  et  al,  1995). 
With  these  constraints,  the  inversion  is  performed  for  the  shear  velocity  of  the  crust  and  upper 
mantle  to  a  depth  of  300  km.  Below  300  km  the  earth  structure  is  fixed,  and  the  inversion  model 
is  required  to  be  continuous  with  the  mantle  structure  at  the  base  of  the  inversion.  In  broad  ocean 
areas,  we  replace  the  Crust  2.0  model  with  models  distinct  for  each  ocean  and  subdivided  by 
ocean  age.  We  also  separate  into  distinct  models  Crust  2.0  models  that  are  geographically 
separated.  So,  for  example,  if  Crust  2.0  has  the  same  model  type  in  North  America  and  in  Asia, 
we  use  the  same  starting  model  for  each,  but  treat  them  as  separate  models  in  the  inversion. 
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2.2  Surface  wave  dispersion  data  set 

During  the  second  year  of  our  project  we  added  three  new  sets  of  dispersion  data  and  improved 
one  already  existing  set  bringing  the  total  number  of  Rayleigh  wave  dispersion  measurements 
with  frequencies  greater  than  or  equal  to  .01  Hz  to  more  than  1,000,000.  One  new  set  comes 
from  Los  Alamos  National  Laboratory  (Yang  et  al,  2002)  consisting  of  more  than  37,000 
dispersion  measurements  (2009  individual  paths)  from  Central  Asia  ranging  between  0.05  Hz 
and  0.23  Hz.  For  frequencies  below  0.03  Hz  these  data  are  much  slower  than  other  data  in  the 
same  region;  consequently  we  have  removed  all  the  Los  Alamos  data  below  0.04  Hz.  Another  set 
comes  from  Huang  et  al  (2003)  and  consists  of  more  than  285,000  data  points  (9730  paths)  from 
China.  These  data  include  unrealistically  fast  paths  crossing  cells  with  water  at  frequencies 
higher  than  about  0.06  Hz.  Consequently  we  have  removed  all  the  Huang  data  points  with  paths 
passing  through  cells  with  water  above  0.04  Hz.  We  have  added  approximately  14,000  our  own 
dispersion  measurements  for  paths  to  IMS  stations  crossing  Eurasia  (Levshin  et  al,  2003).  Pre¬ 
existing  data  coming  from  University  of  Colorado  (Levshin  et  al,  2002)  were  improved  by 
relocating  events  to  the  hypocenters  in  Engdahl  et  al  (1998)  where  possible  as  described  in 
Levshin  et  al.  (2003).  Other  data  already  in  the  data  set  are  described  in  Stevens  et  al  (2001a, b 
and  2002).  Figure  1  shows  the  frequency  distribution  of  group  velocity  and  phase  measurements 
in  our  data  set.  In  addition  to  the  phase  velocity  measurements  shown  in  the  figure,  the  data  set 
also  includes  a  large  number  of  long  period  phase  velocities  derived  from  the  global  phase 
velocity  model  of  Ekstrom  et  al.  (1996). 


Figure  1.  Bar  graph  of  the  number  of  group  velocity  measurements  (left)  and  phase  velocity  measurements  (right) 
in  each  frequency  band  for  all  data  currently  used  in  the  tomographic  inversions. 

2.3  The  inversion  procedure  for  the  3D  earth  model 

The  relationship  between  dispersion  and  the  shear  wave  velocities  of  the  layers  in  the  earth 
model  is  non-linear,  so  the  shear  wave  velocities  are  estimated  by  non-linear  least  squares.  At 
each  step  a  system  of  tomographic  equations  is  formed,  augmented  by  additional  equations  of 
constraint  and  then  solved  by  the  LSQR  algorithm.  The  equations  solved  are 
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where  x  is  the  vector  of  slowness  adjustments  to  the  shear  wave  slownesses  of  layers  in  each  of 
the  572  model  types.  I  is  the  vector  of  slowness  differences  between  predicted  and  observed 
dispersion  measurements,  e  is  the  vector  of  residuals  that  remain  after  inversion  (the  inversion 
minimizes  |e|).  Xo  is  the  vector  of  slownesses  estimated  in  the  last  iteration.  The  elements  of  the 
matrix  A  consist  of  partial  derivatives  of  dispersion  predictions  with  respect  to  shear  wave 
slownesses  in  each  layer.  H  is  a  difference  operator  that  applies  to  vertically  neighboring  layers 
and  has  the  effect  of  constraining  the  vertical  smoothness  of  velocity  profile,  s  is  the  weighting 
of  the  smoothness  constraint  and  can  be  a  diagonal  matrix  (for  variably  weighted  smoothing)  or  a 
scalar  (constant  smoothing).  We  have  implemented  variable  smoothing  so  that  a  different 
smoothing  parameter  can  be  selected  for  each  model  type.  Lateral  smoothing,  which  is  usually 
applied  in  tomography  studies,  is  executed  indirectly  in  our  study  through  selection  of  the  model 
types.  I  is  the  identity  matrix  and  X  weights  the  damping  which  constrains  the  norm  of  the 
difference  between  final  slownesses  and  constraining  model  slownesses  (in  our  work  a  variant  of 
the  Crust  2.0  values).  X  can  be  a  scalar  for  constant  damping,  or  a  diagonal  matrix  for  variable 
damping.  As  for  smoothing,  variable  damping  is  implemented  so  that  a  different  parameter  can 
be  selected  for  each  model  type. 

2.4  Regularization  and  predictions 

Choosing  regularizing  parameters  is  an  essential  part  of  finding  a  model  which  best  predicts 
dispersion,  since  regularization  acts  both  to  control  the  influence  of  data  noise  on  the  estimation 
of  model  parameters  and  to  constrain  parts  of  the  model  that  are  poorly  constrained  by  data.  Too 
much  regularization  will  make  the  model  too  smooth  and  too  little  regularization  will  allow  noise 
to  be  projected  into  the  model,  making  it  unrealistic  and  rough.  In  this  study  the  damping  and 
smoothing  constraints  and  their  associated  weighting  parameters  are  used  to  regularize  the 
solution.  Techniques  for  optimization  of  regularization  parameters  are  not  yet  mature,  especially 
for  large  scale  problems  such  as  this.  The  methods  most  often  described  (e.g.  Hansen,  1998)  are 
the  L-curve,  generalized  cross  validation,  and  discrepancy  principle.  In  the  literature  the  first  two 
methods  are  usually  applied  to  smaller  scale  problems  than  ours  and  with  only  one  regularizing 
parameter,  whereas  we  have  at  least  2  and  possibly  many  more.  The  last  method  mentioned 
requires  a  reliable  independent  estimate  of  data  noise  and  works  by  selecting  the  regularization 
resulting  in  the  residual  based  estimate  of  noise  being  the  same  as  the  independent  estimate.  We 
have  experimented  with  several  of  these  techniques  for  our  inversion  problem,  but  have  not 
found  any  reliable  enough  to  replace  analyst  review  of  the  inversion  results. 

2.5  Techniques  for  evolving  earth  models 

The  selection  of  model  types  and  regularizing  parameters  are  interrelated.  For  example  making  a 
parameterization  finer  (i.e.  adding  new  model  types)  without  increasing  data  coverage  increases 
the  need  for  regularization.  Our  approach  for  finding  a  best  predictive  model  has  been  to  start 
with  a  small  number  of  model  types  and  increase  this  number  gradually  as  new  data  become 
available  or  when  we  detect  systematic  data  misfit.  Once  a  new  parameterization  is  determined, 
the  regularizing  parameters  are  adjusted. 

To  determine  whether  there  are  enough  model  types,  we  carry  out  2D  inversions  of  group 
velocity  residuals  in  at  least  12  narrow  frequency  bands  to  determine  group  velocity  adjustment 
maps  at  1 -degree  resolution.  Each  2D  inversion  of  group  velocity  residuals  is  regularized  with  a 
damping  parameter  which  constrains  the  norm  of  the  group  velocity  adjustment,  and  a  smoothing 
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parameter  which  weights  a  first  difference  operator  that  constrains  lateral  smoothness  of  the 
estimated  group  velocity  (i.e.  v+5v).  We  vary  these  two  parameters  for  each  frequency  band  to 
find  the  smoothest  looking  map  that  still  reduces  data  misfit  reasonably.  The  resulting  maps  are 
examined  to  find  areas  of  similar  adjustment  common  to  most  frequency  bands,  which  are  then 
used  to  delineate  new  model  types. 
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Figure  2.  Shear  velocity  profiles  for  two  model  types:  less  sensitive  to  regularization  on  the  left  and  more  sensitive 
on  the  right.  Ch  is  from  the  Yellow  Sea  and  Max  is  from  North  East  of  Mexico.  The  different  structures 
correspond  to  different  combinations  of  the  scalar  regularization  parameters  s  and  X. 

After  the  model  types  are  selected  the  damping  and  smoothing  parameters  (scalars,  or  diagonal 
matrices)  X  and  s  are  adjusted  to  find  a  reasonable  looking  model  that  still  fits  the  data 
adequately.  Currently  we  are  experimenting  with  spatially  variable  regularization  using  diagonal 
matrices  rather  than  scalars.  We  find  that  the  sensitivity  of  slowness  adjustments  to  scalar 
regularization  parameter  settings  is  not  uniform  and  depends  on  model  type.  In  other  words  when 
comparing  adjustments  for  different  combinations  of  scalar  pairs  (s,  A,),  there  is  much  less 
variation  for  some  model  types  than  others  (see  for  example  Figure  2).  The  strategy  we  are  now 
developing  is  to  relax  regularization  for  those  types  that  are  relatively  insensitive  to 
regularization,  and  to  increase  it  for  those  model  types  that  are  sensitive  to  it. 


A  measure  of  sensitivity  to  regularization  is  the  maximum  of  the  absolute  adjustments  in  S-wave 
slownesses  (max_dev)  for  one  or  more  regularizations.  For  example  the  map  in  Figure  3  shows 
how  this  deviation  varies  geographically  when  (s,  X)=  (50,50).  The  map  shows  that  oceanic 
model  types  are  under  regularized  compared  to  the  others.  Several  regularizations  with  s  and  X 
diagonal  matrices  were  formed  where  (s(i),  A,(i))  =  a*K*max_dev(i),  “i”  indicates  model  type, 
s(i)  and  A,(i)  are  diagonal  matrices  with  dimension  equal  to  the  number  of  adjustable  layers  in  the 
i-th  model,  and  K  is  a  constant.  K  was  varied  to  find  its  best  value  by  comparing  solutions.  The 
vector  a  is  (1,1)  or  (2,1).  After  a  good  value  of  K  is  found  a  small  number  of  model  types  are  still 
too  sensitive,  and  for  them  the  regularization  needed  to  be  increased  further.  Figure  4  shows 
maximum  absolute  deviations  for  the  latest  non-uniform  regularization.  Here  the  absolute 
deviations  are  much  more  similar  among  model  types  than  in  the  uniform  regularization  case. 
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Figure  3.  Map  of  maximum  absolute  slowness  deviations  (units  of  10  for  each  model  t3^e  for  a  uniform 
regularization  of  (s,  X)=  (50,50), 
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Figure  4.  Map  of  maximum  absolute  slowness  deviations  (units  of  10'^)  for  each  model  type  for  most  recent  model 
type  dependent  (non-uniform)  regularization  . 
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2.6  Data  statistics  for  best  current  3D  earth  model 

The  means  and  standard  deviations  of  normalized  group  velocity  residuals,  1-Vo/vp,  where  Voand 
Vp  are  observed  and  model  predicted  group  velocities,  were  calculated  for  narrow  frequency 
bands  and  are  shown  in  Figure  5.  Solid  is  for  our  best  1-degree  model,  and  dashed  is  for  the  best 
5-degree  model  (e.g.  Stevens  and  Adams,  2000)  based  on  Crust  5.1  (Mooney,  et  al.,  1998). 
Figure  5  shows  the  value  of  the  1-degree  model,  especially  for  high  frequencies. 


Figure  5.  Standard  deviations  (o)  and  means  (+)  of  normalized  group  velocity  residuals  are  plotted  against 
frequency  for  Idegree  earth  model  (solid  red)  and  5-degree  earth  model  (dashed  blue). 
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3.  OPTIMIZATION  OF  SURFACE  WAVE  AMPLITUDE  MEASUREMENTS 

Surface  wave  measurements  traditionally  are  accomplished  by  measuring  a  time  domain 
amplitude  at  a  period  near  20  seconds  and  then  calculating  a  surface  wave  magnitude  This 
procedure  is  problematic  at  regional  distances  because  the  surface  wave  is  not  well  dispersed  and 
a  distinct  20-second  arrival  may  not  be  present.  It  is  possible  to  measure  time  domain  amplitudes 
at  higher  frequencies  with  corrections  (e.g.  Marshall  and  Basham,  1972),  however  measurements 
may  be  inaccurate  due  to  differences  in  dispersion  caused  by  differences  in  earth  structure. 
Stevens  and  McLaughlin  (2001)  suggested  as  an  alternative  replacing  time  domain 
measurements  with  a  path  corrected  spectral  magnitude.  The  path  corrected  spectral  magnitude, 
logMo,  is  calculated  by  dividing  the  observed  surface  wave  spectrum  by  the  Green’s  function  for 
an  explosion  of  unit  moment  and  taking  the  logarithm  of  this  ratio,  averaged  over  any  desired 
frequency  band.  The  objectives  of  the  present  study  include  determining  the  optimum  frequency 
band  for  measurement  and  the  best  procedure  for  averaging  the  spectra  over  this  band. 

The  advantages  of  using  logMo  instead  of  the  traditional  surface  wave  magnitude  Ms  are  that 
logMo  is  insensitive  to  dispersion,  independent  of  distance,  works  well  at  regional  distances,  and 
is  inherently  regionalizeable.  Regionalized  path  corrected  spectral  magnitudes  incorporate 
geographic  variations  in  source  excitation  and  attenuation.  Furthermore,  as  discussed  below,  it 
can  in  principle  be  measured  over  different  frequency  bands  to  optimize  the  signal-to-noise  ratio. 
Ms  and  logMo  share  some  limitations:  spectra  from  earthquakes  vary  due  to  source  mechanism 
and  depth,  and  errors  can  occur  if  the  measurement  is  made  in  a  spectral  dip  or  at  high 
frequencies  for  a  deep  event  (Figure  6).  Azimuthal  variations  in  amplitude  caused  by  focal 
mechanism  also  affect  the  amplitudes  of  both  logMo  and  Ms. 


Scalar  Moment  Estimate 


Figure  6.  Path  corrected  spectral  magnitude  for  an  explosion  and  for  earthquakes  calculated  for  several  depths.  The 
path  corrected  explosion  spectrum  is  flat  over  the  entire  frequency  band  (for  perfect  data  and  path 
correction,  while  the  path  corrected  earthquake  spectrum  is  flattened,  but  has  some  variation  due  to 
source  mechanism  and  source  depth. 


9 


The  test  cases  discussed  by  Stevens  and  McLaughlin  (2001)  used  a  frequency  band  of  0.02-0.05  s 
(50-20  s)  to  estimate  the  spectral  magnitudes.  They  estimated  that  on  average,  the  time  domain  and 
spectral  magnitudes  are  related  as  logMo=Ms+11.75.  Most  of  the  waveforms  in  that  work  were 
recorded  at  distances  exceeding  8°.  Due  to  the  relatively  flat  spectra  over  the  0.02-0.05  Hz  band  for 
most  data,  this  choice  worked  quite  well.  The  authors  noted,  however,  that  higher  frequencies 
might  be  required  for  shorter  paths.  An  important  observation  was  that  the  logMo  residuals  are 
independent  of  distance,  despite  the  simple  Q  models  used  in  the  earth  structures  (Figure  7). 


Loq(MO)  vs.  Distance 


Figure?.  Path  corrected  spectral  magnitude  (Log  MO)  residuals  plotted  vs.  distance  (from  Stevens  and 
McLaughlin,  2001).  Log  MO  is  nearly  independent  of  distance. 

In  the  present  work  we  focus  on  the  utility  of  higher  frequencies  in  estimating  spectral 
magnitudes  of  smaller  events,  recorded  at  smaller  distances.  The  purpose  is  to  optimize  the 
spectral  magnitude  estimates,  to  test  their  distance  and  frequency  independence,  and  to  identify 
any  measurement  problems  or  pitfalls.  For  large  amplitude  signals  we  can  expect  the  lower 
frequencies  to  be  better  in  general,  particularly  at  larger  distances  due  to  greater  attenuation  at 
higher  frequencies.  Our  hypothesis  at  the  initiation  of  this  study  was  that  using  higher 
frequencies  for  measuring  spectral  magnitudes  at  shorter  distances  would  optimize  signal  to 
noise  ratio  and  therefore  be  better  for  measuring  surface  waves  at  regional  distances,  however  as 
discussed  below  this  is  only  true  to  a  limited  extent. 

3.1  Surface  wave  measurements  using  events  on  and  near  the  Lop  Nor  test  site 

To  optimize  the  measurement  procedures  and  examine  the  performance  of  the  path  corrected 
spectral  magnitude  at  regional  distances,  we  use  584  spectra  from  76  earthquakes  and  1 1 
explosions  in  the  Lop  Nor  area  (Figure  8).  Additional  spectra  are  available  from  these  events,  but 
only  the  above  spectra  were  deemed  to  be  of  good  quality.  This  means  that  they  passed  the 
dispersion  test  described  by  Stevens  and  McLaughlin  (2001)  and  were  checked  for  certain 
anomalies  such  as  incorrect  instrument  calibrations.  Approximately  11%  of  the  spectra  used  for 
the  logMo  estimates  originate  from  records  at  source-station  distances  of  5°  or  less,  and  another 
11%  at  distances  of  30  or  greater.  Thus  the  bulk  of  the  data  comes  from  intermediate  distances. 
Figure  9  shows  examples  of  explosion  and  earthquake  path  corrected  spectra  from  Lop  Nor  at 
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various  distances.  The  tendency  of  the  explosion  spectra  to  be  relatively  flat  over  more  extended 
frequency  bands  compared  with  earthquake  spectra  is  evident.  This  is  expected  because  1)  the 
spectra  are  corrected  by  an  explosion  Green’s  function  that  flattens  earthquake  spectra 
imperfectly,  and  2)  the  earthquake  spectra  have  frequency  variations  caused  by  source 
mechanism  and  depth. 


Figure  8.  Maps  showing  the  Lop  Nor  area  (rectangle),  stations  (triangles),  and  earthquake  (circles)  and  explosion 
(crosses)  epicenters. 
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Figure  9.  Examples  of  path  corrected  spectra  used  in  this  work:  (a)  Lop  Nor  explosions  recorded  at  distances  of  2®, 
T  and  67“  (left);  (b)  Lop  Nor  earthquakes  recorded  at  distances  of  0.4“,  22“  and  65“  (right).  See  examples 
of  station  logA/o  estimates  in  Table  1.  S/N  is  good  at  all  but  the  highest  frequencies.  Low  frequency  noise 
dominates  over  the  surface  wave  signal  below  about  0.02  Hz. 
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We  calculated  individual  spectral  magnitudes  (i.e.,  several  station  magnitudes  per  event)  over  all 
possible  frequency  bands  between  0.02  Hz  (50  s)  and  0.15  Hz  (~7s),  with  bandwidths  of  0.03 
Hz,  0.04  Hz,  etc.,  up  to  0.13  Hz  for  the  0.02-0.15  Hz  band.  This  procedure  provided  153  bands  to 
examine  from  each  spectrum.  In  the  search  for  the  most  robust  spectral  magnitude  estimate,  four 
different  methods  were  used  as  follows. 

1.  Calculating  a  “simple”  mean  of  the  logarithms  of  all  path  corrected  amplitude 
measurements  made  in  a  given  frequency  band.  This  is  comparable  to  Stevens  and 
McLaughlin’s  (2001)  estimates  in  the  0.02-0.05  Hz  band. 

2.  Iteratively  calculating  a  “robust”  mean,  by  rejecting  outliers  outside  two  standard 
deviations  from  the  mean  calculated  at  each  step.  The  procedure  ends  when  all 
measurements  remain  within  two  standard  deviations  or  when  half  of  the  amplitude 
measurements  in  a  frequency  band  are  rejected,  whichever  occurs  first.  Thus  the  spectral 
magnitude  estimates  are  much  less  affected  adversely  by  the  tendency  of  some  spectra  to 
sharply  vary  in  amplitude  over  some  frequencies,  with  most  outliers  marking 
anomalously  low  amplitudes  (see  Figure  9  above).  Figure  10  compares  the  individual 
(station)  logMo  estimates  from  (1)  and  (2).  Standard  deviations  from  the  robust-mean 
method  are  predictably  lower  than  those  in  the  simple-mean  method,  as  the  insets  in 
Figure  10  show. 

3.  Calculating  logMo  at  the  center  frequency  of  a  least-squares  straight  line  fit  to  the 
spectrum  over  a  given  frequency  band. 


4.  Same  as  (3),  but  the  straight  line  is  “robust”,  minimizing  the  absolute  deviations  of  the 
logarithms  of  the  amplitudes  from  the  line. 


LogMO  /  simple  mean 

Figure  10.  Comparison  of  station  spectral  Figure  11. Histograms  of  standard  deviations  of  the  mean 
magnitudes  calculated  with  two  different  station  logAfo  estimated  in  three  frequency  bands 

methods.  Insets  show  standard  deviations  (shown  on  right),  for  small  and  larger  distances 

as  indicated.  See  text  for  details.  (shown  on  top). 
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The  above  estimates  were  compared  in  order  to  select  the  most  suitable  frequency  bands, 
possibly  varying  with  distance.  Ideally,  the  corrected  spectra  would  be  flat  over  an  extended 
band  of  frequencies.  Flatness  is  particularly  expected  for  explosions,  as  supported  (within  limits) 
by  the  explosion  examples  in  Figure  9  above.  The  magnitude  spectra  estimated  over  any 
reasonable  band  would  be  then  consistent.  In  reality,  truly  flat  spectra  over  extended  frequency 
bands  are  rare,  so  we  need  to  choose  bands  small  enough  not  to  include  too  many  variable 
features  of  the  spectra,  yet  large  enough  not  to  reflect  only  local,  possibly  spurious, 
characteristics. 

In  view  of  the  above,  the  two  main  desirable  properties  of  a  spectrum  over  a  given  frequency 
band  are  small  standard  deviations  and  flatness.  For  this  reason,  in  our  search  for  optimum 
frequency  bands  we  used  two  criteria.  First,  small  standard  deviations  from  (1)  above  represent 
one  measure  of  the  suitability  of  a  frequency  band.  Figure  1 1  indicates  that  for  small  distances, 
the  0.08-0.11  Hz  frequency  band  may  be  preferable  (the  largest  number  of  small  standard 
deviations)  to  either  0.02-0.05  Hz,  or  0.12-0.15  Hz.  Larger  distances  do  not  present  a  clear 
picture,  but  it  is  still  evident  that  relatively  more  small  standard  deviations  are  found  in  the  0.02- 
0.05  Hz  frequency  band,  compared  with  the  higher  frequencies.  We  note  that  at  this  stage  we  do 
not  use  the  standard  deviations  from  (2),  since  they  are  designed  to  greatly  diminish  the  presence 
of  outliers  and  are  thus  not  representative  enough  of  the  quality  of  the  estimates  in  the  different 
frequency  bands.  However,  once  a  suitable  frequency  band  is  chosen,  the  robust  mean  is  the 
most  reliable  estimate  of  logMo.  Spectral  flatness  as  measured  with  the  slopes  of  the  “robust” 
lines  in  (4)  above  provides  a  second  measure  of  the  quality  of  frequency  band;  the  smaller  the 
slope,  the  flatter  the  spectra.  Table  1  shows  examples  of  selected  estimates,  over  one  specific 
frequency  band  out  of  153  (0.05-0.1  Hz),  for  the  explosion  and  earthquake  spectra  shown  in 
Figure  9.  Smaller  slopes  (flatter  spectra)  are  evident  for  explosions  compared  with  earthquakes. 
On  the  other  hand,  increasing  absolute  values  of  slopes  and  standard  deviations  are  seen  for 
earthquakes  with  increasing  distance.  This  is  to  be  expected,  as  the  relatively  high-frequency 
band  in  the  example  is  less  suitable  as  distance  increases. 


Table  1.  Station  LogMO  estimates  in  0.05-0. 10  Hz  from  the  spectra  in  Figure  9. 


Event  Type 

ID.station 

Distance, 

decrees 

Station  logM^ 
(simple) 

Station  logMQ 
(robust) 

Station 

slope/logMo 

mb 

Explosions 

2 1450528.  WMQ 

2.2 

14.31-h0.10 

14.33+0.08 

+1.02/14.31 

gif 

21450535.MAK 

7.1 

15.60+0.15 

15.67+0.05 

+0.41/15.64 

5.8 

21450534.ESDC 

66.8 

14.95+0.22 

15.02+0.09 

+1.62/14.97 

_ _ 

Earthquakes 

21456615.WMQ 

0.4 

13.92+0.14 

14.01+0.05 

+3.64/13.93 

3.2 

21456712.ARU 

122 

14.44+0.25 

14.44+0.24 

-8.70/14.48 

3.8 

21457058.ILAR 

65.3 

15.55+0.27 

15.45-h0.10 

-11.18/15.59 

ESI 
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Figure  12.  Comparison  of  station  spectral  magnitudes  in  different  frequency  bands:  (a)  adjacent  bands,  all  distances; 
(b)  overlapping  bands,  all  distances;  (c)  overlapping  bands,  distances  <5°. 

Next,  we  examined  the  consistency  of  spectral  magnitude  estimates  in  different  frequency  bands. 
Figure  12  shows  examples  of  such  estimates  in  several  frequency  bands  (marked  along  the  plot 
axes).  These  results  indicate  that  although  measurements  are  generally  consistent  when  measured 
in  different  frequency  bands,  some  individual  measurements  do  change  significantly.  Also,  there 
is  a  tendency  for  measurements  to  be  smaller  at  higher  frequencies  (points  lie  slightly  to  the  right 
of  the  lines  in  Figure  12).  These  results  indicate  that  spectral  magnitudes  can  be  measured  in 
different  frequency  bands,  but  with  some  caution  and  attention  to  spectral  shape  variations. 

Finally,  we  examine  which  frequency  bands  perform  best  for  discrimination  between  small 
earthquakes  and  explosions.  That  is,  we  want  to  find  out  if  any  set  of  variable  frequencies  would 
perform  better  in  terms  of  discrimination  than  a  single  frequency  band  applied  at  all  distances. 
Figure  13  shows  logMoimb  plots  using  a  set  of  variable  frequencies  (0.02-0.05  Hz  for  distances 
exceeding  25°,  0.06-0.09  Hz  for  10°  to  250,  and  0.08-0.11  Hz  for  distances  below  10°)  and  the 
0.03-0.07  Hz  frequency  band  for  all  distances.  The  plot  on  the  left,  where  higher  frequencies  are 
used  at  small  distances  (and  hence  for  the  smallest  earthquakes)  apparently  has  a  lower 
discriminating  power  for  small  events  than  when  0.03-0.07  Hz  magnitudes  are  used.  The  reason 
is  that  the  spectral  magnitudes  of  smaller  events  (logMo  14  to  15;  i.e.,  Ms  2.2  to  3.2),  recorded 
predominantly  at  regional  distances,  are  generally  smaller  than  the  estimates  at  lower 
frequencies.  We  examined  the  logMoimb  ratio  for  a  number  of  frequency  bands  and  established 
that  the  0.03-0.07  Hz  interval  performs  best  in  discriminating  between  earthquakes  and 
explosions  for  the  Lop  Nor  data  set.  The  performance  of  the  time  domain  Ms:mb  discriminant  for 
these  events  (not  shown),  is  similar  to  that  of  the  variable  frequency  measurement  (Figure  13, 
left),  although  the  comparison  is  complicated  by  the  fact  that  there  is  not  a  standard  procedure 
for  measuring  time  domain  Ms  in  the  regional  distance  range. 
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Figure  13.LogMo:mb  plots  showing  event  spectral  magnitudes  for  earthquakes  (O)  and  explosions  (X)  in  Lop  Nor. 

Station  spectral  magnitudes  were  calculated  using  frequencies  increasing  with  distance  (left)  and  the 
0.03-0.07  Hz  frequency  band  for  all  distances  (right).  Bold  lines  mark  the  empirical  discrimination 
relationship  of  Stevens  and  McLaughlin  (2001). 

3.2  Path  corrected  noise  estimates 

While  a  path  corrected  spectral  magnitude  can  be  measured  over  any  frequency  band,  it  is 
subject  to  the  following  constraints: 

1.  Earthquake  spectra  decrease  at  high  frequencies,  depending  on  depth  (Figure  6),  so  the 
high  end  of  the  frequency  band  should  be  low  enough  that  discrimination  is  not  degraded 
by  this  effect  (Figure  13). 

2.  At  high  frequencies,  attenuation  is  higher  and  the  dispersion  more  variable,  so  the  path 
correction  is  likely  to  be  better  and  the  signal  may  be  higher  at  lower  frequencies. 

3.  Noise  increases  at  low  frequency,  so  the  low  end  of  the  frequency  band  should  be  at  a 
frequency  high  enough  to  be  above  the  noise  level. 

In  order  to  quantify  these  effects  better,  we  measured  some  “path  corrected  noise  spectra”  for  the 
Lop  Nor  data  set.  These  are  simply  noise  spectra  measured  at  the  station  that  have  been  divided  by 
an  explosion  Green’s  function  in  the  same  manner  as  would  be  done  for  a  signal.  Since  the  signal 
spectra  are  approximately  flat  over  most  of  the  frequency  band,  the  path  corrected  noise  spectra  are 
a  measure  of  the  minimum  path  corrected  signal  that  could  be  measured  at  each  station. 

Figure  14  shows  the  noise  measurements  for  two  stations;  WMQ,  located  an  average  of  2 
degrees  from  each  event,  and  HIA,  located  an  average  of  27  degrees  from  each  event.  At  the 
lowest  frequencies,  the  noise  levels  decrease  with  increasing  frequency,  reaching  a  minimum  at 
about  0.04  Hz,  then  increase  to  a  peak  at  about  0.06  Hz,  which  corresponds  to  the  frequency  of 
the  primary  microseism  peak.  The  noise  levels  then  decline  to  about  0.1  Hz,  and  then  increase. 
As  Figure  14  shows,  there  is  also  substantial  variability  in  the  noise  level,  so  that  although  the 
average  noise  level  suggests  that  the  minimum  spectral  magnitude  that  could  be  measured  is 
about  13.5  at  WMQ  and  14.5  at  HIA,  it  may  be  substantially  lower  or  higher.  These  results 
suggest  that  in  general  it  is  best  to  measure  surface  wave  spectra  at  frequencies  above  0.03  Hz 


15 


and  below  0.1  Hz.  Although  noise  levels  remain  fairly  low  in  the  0.1 -0.2  Hz  frequency  band  at 
the  closest  stations,  factors  #1  and  #2  mentioned  above  make  measurement  in  this  band  risky. 
Based  on  the  empirical  tests  described  earlier,  we  recommend  measuring  surface  waves  at 
frequencies  between  0.03  and  0.08  Hz. 


Path  Corrected  Noise  Spectra  at  HI  A 


Figure  14.  Average  path  corrected  noise  measurement  and  ±one  standard  deviation  curves  for  13  time  segments  at 
WMQ  (bottom)  and  for  54  time  segments  at  HIA  (top).  The  average  distance  to  WMQ  is  2  degrees  and 
the  average  distance  to  HIA  is  27  degrees. 
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4.  IMPROVED  AZIMUTH  ESTIMATION 


To  monitor  nuclear  testing,  the  International  Data  Center  (IDC)  automatically  processes  surface 
waves  recorded  at  individual  International  Monitoring  System  (IMS)  stations,  attempting  to 
identify  them  and  measure  amplitude  and  Mj.  Three-component  azimuth  estimates  derived  using 
the  current,  spectral  method  are  recorded  but  are  too  inaccurate  to  use  for  identification.  That 
method  is  based  on  an  algorithm  originally  proposed  (Smart,  1978)  as  part  of  a  surface  wave 
detector,  rather  than  as  a  means  of  estimating  the  azimuth.  Selby  (2000)  suggested  that  another 
detection  technique  (Chael,  1997),  based  on  correlation  of  the  vertical  and  Hilbert  transformed 
radial  component  Rayleigh  wave  records,  could  be  used  for  azimuth  estimation  and  would 
improve  IDC  processing. 

To  compare  the  performance  of  the  algorithms,  we  incorporated  the  Chael/Selby  (CS)  algorithm 
into  the  automatic  surface  wave  processing  software  that  utilizes  the  current  IDC  method  (CM). 
In  the  standard  processing  of  surface  waves,  there  is  a  detection  if  the  surface  waves  pass  a 
dispersion  test  (Stevens  and  McLaughlin,  2001).  In  this  study,  both  algorithms  are  tested  on 
2,599  records  that  passed  the  dispersion  test  and  2,363  generally  poorer  S/N  records  that  failed, 
all  from  events  that  had  at  least  four  surface  wave  detections.  We  assess  the  performance  of  each 
algorithm  for  different  passbands  and  group  velocity  windows. 


4.1  Differences  between  the  algorithms 

Both  the  CM  and  CS  algorithms  find  the  backazimuth  that  best  matches  the  expectation  that  the 
horizontal  and  vertical  components  of  the  Rayleigh  wave  records  are  similar,  but  90°  out  of 
phase.  The  biggest  difference  between  algorithms  is  that  the  CM  uses  Love  as  well  as  Rayleigh 
waves  to  estimate  the  backazimuth. 


The  use  of  both  Love  and  Rayleigh  waves  by  the  CM  versus  only  Rayleigh  waves  by  the  CS 
algorithm  makes  it  impossible  to  compare  the  algorithms  using  the  same  time  window  while 
optimizing  the  performance  of  both.  The  CS  algorithm  should  perform  optimally  when  the  time 
window  encompasses  just  the  Rayleigh  wave,  as  a  longer  window  will  only  add  noise. 
Antithetically,  the  CM  should  perform  optimally  with  a  longer  time  window  that  encompasses 
both  the  Rayleigh  wave  and  the  higher  group  velocity  Love  wave. 


4.1.1  The  Current  Method 

The  CM  finds  values  for  rn  and  L,  the  complex  Fourier  coefficients  of  the  Rayleigh  and  Love 
wave  vertical  component  displacements,  ((),  the  backazimuth,  and  En,  the  Rayleigh  wave 
ellipticity,  that  minimize  the  squared  distance  between  the  data  and  model  (Smart,  1978;  1981). 
The  ellipticity  may  be  assumed  known,  which  reduces  by  one  the  degrees  of  freedom.  Subscripts 
refer  to  frequency.  The  function  minimized  is  written 


\Zn-r„\  + 


yn-(i^nr„  oos{(^)-l„  sin((Z>))|"  +\x„  -{ie„rj  +  l„a)f 


(1) 


where  Xn,  yn,  and  Zn  are  the  complex  Fourier  coefficients  of  the  east,  north,  and  vertical 
components  of  the  seismic  record.  Two  important  elements  of  this  function  are  that  the  radial 
component  of  the  Rayleigh  wave  equals  iEnrn,  that  is,  the  vertical  component  advanced  by  90° 
and  scaled  by  the  ellipticity,  and  that  the  Love  wave  is  independent  of  the  Rayleigh  wave. 
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4.1.2  The  Chael/Selby  algorithm  and  its  implementation 

The  CS  algorithm  (Chael,  1997;  Selby,  2000)  finds  the  backazimuth  for  which  the  correlation 
coefficient  of  the  vertical  and  Hilbert  transformed  radial  is  a  maximum  (Equation  2). 


where, 

N 

Sjk=^XjiT)x^{T).  (3) 

r=l 

The  implementation  uses  the  correlation  of  the  Hilbert  transformed  vertical  with  the  radial,  for 
one  degree  increments  of  backazimuth.  While  Selby  uses  the  maximum  positive  correlation  as 
the  backazimuth,  Chael  (1997)  uses  the  central  azimuth,  determined  by  the  circular  mean  (e.g. 
Fisher,  1993).  We  test  the  algorithm  using  both  measures. 

For  synthetic  Rayleigh  waves  with  azimuthally  evenly  distributed  random  noise,  the  circular 
mean  provides  more  accurate  estimates  than  the  maximum  correlation.  The  lower  the  S/N,  the 
greater  the  advantage  of  the  circular  mean. 

For  noise  free  data  (i.e.  synthetics),  equation  2  will  return  a  negative  or  positive  constant,  as  the 
cross-correlation  of  the  Hilbert  transformed  vertical  and  radial  and  the  autocorrelation  of  the 
radials  change  in  S3mch  as  the  algorithm  steps  through  backazimuths.  Our  implementation  of  the 
algorithm  therefore  provides  two  methods  of  estimating  the  similarity  of  the  components.  One 
uses  the  correlation  coefficient  (equation  2).  The  second  method  is  intended  to  avoid  the 
problem  described  above  when  data  are  noise  free  by  normalizing  by  alone  (equation  4). 


(4) 


4.2  Performance  of  the  algorithms  in  4  frequency  hands. 

The  first  set  of  tests  used  a  2.5  to  5.0  km/sec  group  velocity  window,  which  should  encompass 
both  Love  and  Rayleigh  waves,  and  compared  backazimuth  estimates  in  3  relatively  narrow, 
overlapping  frequency  bands,  0.02-0.04  Hz,  0.03-0.05,  0.04-0.06  Hz,  plus  one  frequency  band 
covering  the  entire  spectrum  from  0.02  to  0.06  Hz. 

Figure  15  shows  histograms  of  azimuth  residuals  for  the  CM  and  one  implementation  of  the  CS 
algorithm  for  the  middle  frequency  band.  Table  2  presents  the  errors  associated  with  each 
method.  The  CS  algorithm  performs  better  in  two  ways.  First,  the  CM  applied  as  it  is  currently 
used  in  the  automatic  processing,  has  sign  errors  for  a  significant  number  of  events.  Second,  the 
histograms  reveal  a  much  larger  number  of  outliers  for  the  CM  compared  to  the  CS  algorithm. 

Surprisingly,  the  best  implementation  of  the  CS  algorithm  is  the  combination  of  the  peak 
position  of  the  correlation  and  normalization  of  the  correlation  by  .  Why  the  maximum 
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correlation  provides  a  more  accurate  estimate  of  backazimuth  than  the  circular  mean  bears 
further  investigation.  One  possibility  is  that  the  real  noise  is  azimuthally  unevenly  distributed  and 
so  biases  the  circular  mean  measurement  toward  the  direction  of  the  predominant  noise  source. 
Such  skewed  noise  would  not  so  strongly  affect  the  maximum  correlation  position. 


Current  I  DC  Method 


Azimuth  Residual  (degrees) 


Selby-Chael  Algorithm 


-100  0  100 

Azimuth  Residual  (degrees) 


Figure  15.  Azimuth  residuals  for  0.03-0.05  Hz  using  the  current  method  (left)  and  the  Chael/  Selby  algorithm  using 
the  normalization  of  equation  2  and  peak  correlation  (right). 

The  scaled  median  absolute  deviation  (SMAD),  the  one-norm  measure  of  the  central  tendency, 
provides  a  more  accurate  measure  of  the  spread  of  data  about  the  central  value  for  such  heavy¬ 
tailed  data  as  these  than  does  the  standard  deviation  (STD).  Errors  for  all  of  the  tests  performed 
badly  fail  the  Kolmogorov-Smimov  test  for  normality  of  a  distribution.  That  the  STDs  are  so 
much  larger  than  the  SMADs  indicates  that,  while  most  estimates  are  reasonably  accurate,  there 
are  a  large  number  of  outliers. 

The  CM  performs  best  at  higher  frequency,  while  the  CS  algorithm  performs  best  at  lower 
frequency.  The  CM  performs  almost  as  well  as  the  CS  algorithm  at  the  highest  frequencies 
tested.  This  is  likely  because  the  Love  waves  are  relatively  larger  at  longer  periods  than  the 
Rayleigh  waves,  and  contrary  to  expectations,  the  CM  performs  poorly  when  the  Love  waves  are 
very  large.  This  is  probably  due  to  the  180  degree  ambiguity  in  Love  wave  polarization,  which 
causes  large  Love  wave  amplitudes  to  increase  the  likelihood  of  a  180  degree  error  in  azimuth. 

Larger  azimuth  residuals  at  higher  frequency  for  the  CS  method  may  be  due  to  lower  S/N. 
Alternately,  greater  lateral  heterogeneity  in  the  shallower  Earth  layers  could  lead  to  more 
variation  in  the  propagation  direction  of  higher  frequency  surface  waves.  The  estimates  may  in 
that  case  accurately  reflect  the  effect  of  lateral  heterogeneities  on  higher  frequency  waves. 

For  the  CS  algorithm,  the  backazimuth  residuals  for  the  broader  frequency  band  are  comparable 
to  those  in  the  0.03-0.05  Hz  passband,  larger  than  those  of  the  lowest  frequency  passband,  and 
smaller  than  those  of  the  highest  frequency  passband.  The  parts  of  the  frequency  band  where  the 
backazimuth  either  varies  or  is  less  well  resolved  appear  to  diminish  resolution  of  the  broadband 
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estimate  compared  to  narrower  passband  estimates,  so  a  broad  frequency  band  does  not  provide 
any  advantage.  For  the  CM,  the  performance  is  actually  poorer  for  the  broadband  data  than  for 
each  of  the  narrower  passbands  used.  The  rest  of  the  analyses  are  performed  for  just  the  three 
narrow  frequency  bands. 

Table  2.  Scaled  median  absolute  deviations  and  standard  deviations  (in  parentheses)  of  the  azimuth  residuals  in  4 
frequency  bands  from  the  algorithms  tested.  Four  variations  of  the  implementation  of  the  CS  algorithm  are 
tested.  CS,  and  CS2  use  the  peak  amplitude  of  the  correlation  while  CS3  and  CS4  use  the  circular  mean.  CS, 

and  CS3  are  normalized  by  yjS^Sj^  as  in  equation  2.  CS2  and  CS4  are  normalized  by  as  in  equation  4. 


.02-.04  Hz 

.03-.05  Hz 

.04-.06  Hz 

.02-.06  Hz 

CM 

31 (85) 

22  (71) 

22  (64) 

20  (68) 

CSl 

13(41) 

16  (47) 

19  (53) 

16  (45) 

CS2 

17  (40) 

19  (46) 

22  (52) 

20  (45) 

CS3 

15  (39) 

17(46) 

20  (52) 

17(45) 

CS4 

17  (40) 

19  (46) 

22  (53) 

P20  (45) 

4.3  The  effect  of  signal  strength  (event  size)  on  backazimuth  estimates. 

Figure  16  shows  the  backazimuth  residuals  for  each  of  the  algorithms  vs.  event  size,  which 
serves  as  a  proxy  for  the  signal-to-noise  ratio.  The  CM  is  almost  comparable  in  performance  to 
the  CS  algorithm  for  large  events,  except  for  a  large  number  of  outliers  at  180  degrees.  It  begins 
to  fail  badly  for  smaller  events,  especially  at  the  lowest  frequencies.  Table  3  show  results 
comparable  to  Table  2,  but  for  events  of  Ms  5.0  and  above. 


Table  3.  Same  as  Table  2,  but  for  335  records  of  events  with  Ms  >  5.0. 


.02-.04  Hz 

.03-.05  Hz 

.04-.06  Hz 

CM 

14.5  (68.6) 

14.4  (51.9) 

12.9  (47.3) 

SC, 

8.7  (30.5) 

9.6  (29.8) 

10.6  (35.2) 

SC2 

10.8  (29.5) 

13.0  (28.4) 

13.4  (34.5) 

SC3 

9.4  (29.7) 

11.7  (28.9) 

12.2  (34.9) 

SC4 

10.7  (29.5) 

12.8  (28.4) 

13.2  (34.5) 
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0.02  -  0.04  Hz 


0.04  -  0.06  Hz 


Figure  16.  Median  ±2  SMAD  confidence  intervals  of  azimuth  residuals  binned  by  Ms  values,  for  the  CM  (dotted) 
and  the  best  performing  implementation  of  the  CS  algorithm  (solid).  Each  bin  has  approximately  250 
azimuth  residuals.  Results  in  the  20-33  second  period  passband  are  intermediate  between  those  shown. 
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4.4  Variable  group  velocity 

The  work  reported  above  used  2.5  to  5.0  km/sec  group  velocity  windows  that  we  expect  to  favor 
the  CM,  as  they  are  most  likely  to  include  both  Love  and  Rayleigh  waves.  Next  we  use  the  group 
velocities  predicted  by  a  1°  global  surface  wave  model  (Stevens,  et  al  2002)  to  select  windows 
intended  to  more  narrowly  bracket  the  Rayleigh  waves  and  so  minimize  the  effect  of  noise 
outside  those  windows  (Figure  17).  The  minimum  group  velocity  is  that  used  in  the  dispersion 
test  for  surface  waves  at  the  high  frequency  comer  of  the  passband.  Similarly,  the  maximum 
group  velocity  is  determined  by  the  low  frequency  comer. 

Table  4  is  similar  to  Table  2,  but  for  the  variable  length  group  velocity  windows.  The  performance 
of  the  CM  is  slightly  worse  for  the  shorter  windows  at  the  highest  frequency  passband,  as  expected, 
and  there  is  a  small  improvement  in  performance  for  the  best  performing  algorithm,  CS  using  the 

peak  of  the  correlation  with  normalization  by  yJS^S-  ,  compared  to  the  longer  group  velocity 

windows.  Overall,  the  effect  of  narrowly  isolating  the  Rayleigh  waves  is  quite  small,  which  is 
encouraging  regarding  the  prospects  of  the  CS  algorithm  for  routine  detection. 


Figure  17.  Surface  wave  records  in  the  3  passbands  tested  at  the  station  ARCES  for  an  Ms  5.0  event  at  6035  km. 

The  entire  trace  is  the  5.0  to  2.5  km/s  group  velocity  window,  while  the  outlined  segments  are  the  shorter 
group  velocity  windows  used  to  isolate  the  Rayleigh  waves.  The  isolation  works  best  at  the  highest 
frequency  passband,  where  the  large  Love  wave  (top  trace)  is  outside  the  narrower  window. 

Table  4.  Same  as  Table  2,  but  for  group  velocity  windows  designed  to  isolate  the  Rayleigh  waves.  The  change 

from  the  2.5  to  5  km/s  group  velocity  windows  is  given  in  parentheses  (negative  change  is 
improvement). 


.02-.04  Hz 

.03-.05  Hz 

.04-.06  Hz 

CM 

31.4  (+0.7) 

21.2  (-0.4) 

23.9  (+1.5) 

SC, 

12.7  (-0.5) 

14.9  (-0.8) 

18.1  (-0.5) 

SC2 

17.0  (+0.5) 

19.0  (+0.3) 

21.5  (-0.3) 

SC3 

14.7  (-0.1) 

16.8  (-0.2) 

19.6  (-0.4) 

SC4 

16.9  (+0.4) 

19.0  (+0.2) 

21.6  (-0.2) 
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4.5  Predicting  accuracy  of  estimates  using  the  crosscorrelation  value 

The  crosscorrelation  of  the  Hilbert  transformed  vertical  with  the  radial  predicts  the  accuracy  of 
the  azimuth  estimate.  Figure  18  shows  azimuth  estimate  residuals  vs.  the  cross-correlation  for  the 
best  implementation  of  the  CS  algorithm  applied  to  data  that  pass  the  dispersion  test. 

For  0.02  to  0.04  Hz,  60%  of  the  data  have  a  cross-correlation  >  0.8.  The  SMAD  of  the  error  of 
those  data  is  8.5°,  vs.  12.7°  for  all  the  data.  For  0.03  to  0.05  Hz,  58%  of  the  data  have  a 
crosscorrelation  value  >  0.8,  and  those  data  have  a  smad  of  9.3°  vs  14.9°  for  all  the  data.  For  0.04 
to  0.06  Hz,  54%  of  the  data  have  a  crosscorrelation  value  >  0.8,  and  those  data  have  a  smad  of 
11.1°  vs  18.1°  for  all  the  data. 

This  ability  to  predict  backazimuth  estimate  accuracy  can  aid  association  with  known  events.  In 
particular,  the  strict  dispersion  criteria  for  detection  at  the  IDC  could  be  relaxed  in  cases  where 
the  backazimuth  is  consistent  with  the  theoretical  backazimuth  and  the  crosscorrelation  value  is 
high.  Well  constrained  uncertainties  are  also  important  for  understanding  the  accuracy  of 
inferences  made  from  measurements  of  the  surface  waves,  for  example,  tomography  based  on 
surface  wave  polarization  (e.g.  Yoshizawa,  et  al.,  1999). 


Figure  18.  Median  azimuth  residuals  ±2  SMAD  confidence  intervals  vs.  the  median  cross-correlation  of  the  radial 
and  Hilbert  transformed  vertical  Rayleigh  waves,  for  each  bin  of  171  measurements,  for  the  0.03-0.05  Hz 
passband.  Results  are  similar  for  the  other  passbands. 

4.6  Events  for  which  surface  waves  are  not  detected  . 

Figure  19  shows  the  results  of  each  method  of  polarization  analysis  applied  to  the  poorer  S/N 
records  that  did  not  pass  the  dispersion  test.  The  CS  method  (implemented  using  maximum 
cross-correlation)  extracts  accurate  backazimuth  estimates  for  many  of  these  data.  The  more 
accurate  estimates  can  be  identified  by  their  cross-correlation  value. 


23 


Figure  19.  Azimuth  residuals  for  the  0.02  -  0.04  Hz  passband  for  the  CS  algorithm  (top  row)  and  the  CM  (bottom 
row).  Histograms  of  azimuth  residuals  for  all  the  data  (left  column)  and  for  data  with  correlation  >  0.8  for 
CS  (upper  middle)  or  F-statistic  >  30  for  the  CM  (lower  middle).  The  central  plots  use  -14.5%  of  the 
data  in  each  case.  The  upper  right  plot  shows  the  median  azimuth  residual  ±2  SMAD  confidence 
intervals  vs.  the  correlation,  as  in  Figure  18,  for  the  CS  algorithm.  The  bottom  right  shows  the  median 
azimuth  residual  ±2  SMAD  confidence  vs.  the  F  statistic  of  the  CM. 

The  SMAD  of  the  azimuth  residuals  for  data  not  passing  the  dispersion  test  is  58°,  vs.  108°  for 
the  CM.  The  upper  central  plot  shows  histograms  just  for  data  with  a  correlation  >  0.8,  which 
comprise  14.7%  of  the  data.  Their  azimuth  residual  has  a  SMAD  of  14.8°.  The  F-statistic  of  the 
CM  does  not  provide  similar  predictive  capabilities.  In  fact,  for  the  largest  F-statistic  values,  the 
CM’s  error  increases.  This  is  because  the  F-statistic  can  be  very  large  when  the  signal  is 
dominated  by  very  large  Love  waves,  but  such  records  often  have  180°  (i.e.  sign)  errors.  Even 
disregarding  the  sign  errors,  the  F-statistic  is  a  much  poorer  predictor  of  the  accuracy  of 
polarization  estimates,  and  the  resolution  of  the  current  method  is  much  poorer  than  that  of  the 
CS  algorithm. 

4.7  Conclusions  of  Azimuth  Study 

The  correlation  of  the  radial  and  Hilbert  transformed  vertical  seismic  records  provides  accurate 
estimates  of  Rayleigh  wave  polarization  and  provides  significant  improvement  over  the  current 
method  employed  at  the  IDC.  Further,  the  value  of  the  cross-correlation  provides  a  reliable 
estimate  of  the  accuracy  of  the  polarization. 

The  CS  method  could  also  improve  detection  of  Rayleigh  waves,  through  better  association  with 
known  events.  Specifically,  current  strict  IDC  requirements  on  detection  in  multiple  passbands 
could  be  relaxed  when  there  is  a  backazimuth  estimate  that  has  a  high  cross-correlation  value  and 
is  consistent  with  the  theoretical  backazimuth. 
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5.  CONCLUSIONS  AND  RECOMMENDATIONS 


Improvements  to  surface  wave  dispersion  models  are  being  accomplished  by:  1)  continuing 
addition  of  new  data  with  good  quality  control  and  particular  attention  to  regions  with  gaps  in 
data  coverage;  2)  removal  of  poor  quality  data  in  the  data  set  (poor  quality  data  becomes  more 
apparent  as  the  data  set  increases  and  the  model  quality  improves);  3)  addition  of  model  types 
where  the  data  requires  them;  4)  improvement  in  constraints  on  sediments  and  Moho  thickness; 
5)  improvement  to  the  regularization  techniques,  which  can  now  be  defined  on  a  model  by  model 
basis,  allowing  improved  data  fit  while  achieving  realistic  earth  models. 

Improved  methods  for  surface  wave  measurement  are  being  implemented  and  tested.  Surface 
wave  spectra  are  derived  from  phase-matched  filtered  data,  and  the  phase-matched  filters  are 
derived  from  the  regionalized  dispersion  models.  Path  corrected  spectral  magnitudes  are  derived 
by  dividing  the  observed  spectra  by  an  explosion  Green’s  function,  where  the  Green’s  functions 
are  calculated  from  the  global  earth  models.  Thus  we  use  the  earth  and  dispersion  models  to 
optimize  spectral  measurements  and  regionalize  surface  wave  excitation  and  attenuation.  In  this 
paper,  we  have  described  a  detailed  study  of  procedures  for  optimizing  measurement  of  path 
corrected  spectral  magnitudes.  A  significant  advantage  of  logMo  over  Ms  is  that  it  can  be 
measured  at  any  distance  range  without  the  anomalies  caused  by  variations  in  dispersion  that 
affect  Ms.  In  principle,  logMo  can  be  measured  over  any  frequency  band  and  optimized  by 
choosing  the  band  with  maximum  S/N.  However,  we  found  that  logMo  for  earthquakes  is 
frequently  significantly  lower  at  higher  frequencies,  which  degrades  discrimination,  and  that 
furthermore  the  S/N  for  lower  frequencies  is  good  even  for  very  short  distances.  We  therefore 
recommend  that  surface  wave  measurements  be  made  at  lower  frequencies  even  at  short 
distances.  We  are  in  the  process  of  determining  the  optimum  frequency  band  for  measurement, 
and  our  current  recommendation  is  to  use  a  frequency  band  of  0.03-0.07  Hz  consistently  for  all 
data.  We  are  continuing  to  evaluate  this  recommendation  for  a  larger  data  set  with  more  types  of 
earth  structure. 
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6.  DATA  DELIVERABLE 


Two  data  deliverables  are  provided  together  with  this  annual  report.  They  are: 

1.  The  current  set  of  earth  models,  dispersion  curves  and  other  derived  data.  These  are 
contained  in  compressed  tar  file  “LP_2003_Oct.tar.Z”.  Information  on  data  format  is 
included  in  the  delivery. 

2.  The  maxpmf  program  compiled  for  Sun  Solaris,  and  the  maxpmf  man  page.  These  are 
contained  in  compressed  tar  file  “maxpmf_5.2.tar.Z” 

These  data  deliverables  may  be  obtained  from  the  contracting  office  or  from  the  authors. 
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